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Abstract. The need to compute small con-eigenvalues and the associated con-eigenvectors of 
positive-definite Cauchy matrices naturally arises when constructing rational approximations with 
a (near) optimally small L°° error. Specifically, given a rational function with n poles in the unit 
5h , disk, a rational approximation with m <^n poles in the unit disk may be obtained from the mth 

con-eigenvector of an n X n Cauchy matrix, where the associated con-eigenvalue Am > gives 
the approximation error in the L°° norm. Unfortunately, standard algorithms do not accurately 
compute small con-eigenvalues (and the associated con-eigenvectors) and, in particular, yield few 
or no correct digits for con-eigenvalues smaller than the machine roundoff. We develop a fast 
and accurate algorithm for computing con-eigenvalues and con-eigenvectors of positive-definite 
Cauchy matrices, yielding even the tiniest con-eigenvalues with high relative accuracy. The algo- 
rithm computes the mth con-eigenvalue in O (m?n) operations and, since the con-eigenvalues of 
positive-definite Cauchy matrices decay exponentially fast, we obtain (near) optimal rational ap- 
proximations in O (log5~^)'^^ operations, where 5 is the approximation error in the L°° norm. 
We derive error bounds demonstrating high relative accuracy of the computed con-eigenvalues 
and the high accuracy of the unit con-eigenvectors. We also provide examples of using the al- 
gorithm to compute (near) optimal rational approximations of functions with singularities and 
sharp transitions, where approximation errors close to machine precision are obtained. Finally, we 
present numerical tests on random (complex-valued) Cauchy matrices to show that the algorithm 
computes all the con-eigenvalues and con-eigenvectors with nearly full precision. 



1. Introduction 

We present an algorithm for computing with high relative accuracy the con-eigenvalue decompo- 
sition of positive-definite Cauchy matrices, 



(1.1) CUm = X,nUm, Cij = ' «, j = 1, . . . , n, 

where 7^ and are complex numbers and |7i| < l.The con-eigenvalue A„j is only defined up to an 
arbitrary phase, which we choose so that A,„ > 0. Although the con-eigenvalue decomposition (see 
e.g. |30| ) is less well-known than the eigenvalue decomposition or the singular value decomposition, 
it arises naturally in constructing optimal approximations using exponentials or rational functions 
[11 [21 [21 [131 US [SI [7] . For example, for a real- valued rational function f{z), 



(1.2) /w = E^+Et^+"o, 



1=1 i=i 



we may construct a rational approximation g{z) with m poles and with an error, 



max 

a;e[0,l] 



by solving the con-eigenvalue problem (|l.ip (see Section 12.11 for more detail) . Ordering the con- 
eigenvalues, Ai > . . . > A„ > 0, the number of poles m of the approxiniant g{z) corresponds to the 
index of the con-eigenvalue Am and leads to a near optimal approximation in the L°°-norm with 



This research was partially supported by NSF grant DMS-100995 and DOE/ORNL grant 4000038129. 



1 



CON-EIGENVALUE ALGORITHM FOR OPTIMAL RATIONAL APPROXIMATIONS 



2 



the error close to Am- The form (|1.2p ensures that / (e^'^") is real- valued and periodic; complex- 
valued functions may also be handled using this form by splitting the real and imaginary parts and 
performing additional reductions (see [7J). 

Current algorithms compute an approximate con-eigenvalue A,„ with an error no better than 
Am — Am / |Ai| — O (e), and an approximate unit con-eigenvector with an error no better than 

l!"m -11^112 = O (e) /absgap,„, absgap„ = min |Am - Ap| / |Ai| , 

where e denotes the machine roundoff. This implies that a computed con-eigenvalue smaller than 
|Ai| e may have few or no correct digits. Hence, in order to obtain a rational approximation with 
accuracy Am ^ 10"'', we may be forced to use at least quadruple precision. Since quadruple precision 
is typically not supported by the hardware, it slows down the computation by an unpleasant factor 
(between 30 and 100). Another undesirable feature of current algorithms to solve (jl.ip is the O (n^) 
complexity for finding the m <^ n poles of g{z), where n is the original number of poles of f{z). 

Although the construction of optimal rational approximations in the L°°-norm has a long history 
(starting with the seminal papers [3 [H [3]), the difficulties mentioned above limit practical appli- 
cations of such approximations to situations where the problem size is relatively small and a low 
accuracy is acceptable. In this regard, we view our results as a stepping stone toward a wider use 
of optimal L°°-approximations in numerical analysis (see \27\]. 

We develop a fast and accurate algorithm for con-eigenvalue/con-eigenvector computations of 
positive-definite Cauchy matrices that addresses both of the difficulties mentioned above. Our 
algorithm computes the mth con-eigenvalue/con-eigenvector in O (rn^rij operations (see Section[5]). 
Since the con-eigenvalues of positive definite Cauchy matrices decay exponentially fast, for a given 
desired accuracy ||/ (e^'^"') — g (e^'^™) ||oo ~ the number of poles m in the approximant g{z) is 

O (logS^^y Therefore, the complexity of our algorithm is O (log (5^^)^^ , i.e., it is essentially 
linear in the number of original poles n and, thus, is mostly controlled by the number of poles of 
the final optimal approximation. 

The con-eigenvalue algorithm achieves high relative accuracy, i.e., the computed con-eigenvalue 



Am satisfies 



Am, ^ A^ 



/ I Am I = O {e), and the computed unit con-eigenvector satisfies 



ljum - u^i\\2 O (e) /relgap„, relgap,„ = min |Am - A;| / (A; + Am) : 



(see Theorems [6] and [7] for the exact statement) . In contrast to the usual perturbation theory for 
general matrices, we show that small perturbations of the poles 7m and residues am (determining 
the Cauchy matrix C = C{a,j) in (|l.ip ) lead to correspondingly small perturbations in the con- 
eigenvalues and con-eigenvectors, as long as the poles are well separated in a relative sense and are 
not too close to the unit circle. 

In many applications, the function / (e^^") has sharp transitions, so that the poles are clustered 
close to the unit circle and each other. In such cases, it is natural to maintain the poles of / (z) 
in the form = cxp(— r^), where TZe (tj) > and < Xm{Tj) < 2tt, so that TZe{Tj) are well- 
separated in a relative sense. The reduction algorithm produces new poles of the same form, where 
even the smallest exponents are computed with high relative accuracy. This allows us to develop a 
numerical calculus that includes functions with singularities and sharp transitions. We address this 
issue further in Section |3l 

Our approach is inspired by papers [20l [23l [HI [TH [28] , which develop algorithms and theory for 
highly accurate SVDs of certain structured matrices. Generally speaking, high relative accuracy is 
achieved when it is possible to avoid catastrophic cancellation resulting from subtracting two close 
floating point numbers (when the outcome of such cancellation is significant relative to the final 
result). We refer to [16] for a comprehensive analysis of when efficient and accurate algorithms 
are possible using floating point arithmetic. Classes of matrices for which highly accurate SVD or 
eigenvalue algorithms exist include bi-diagonal matrices |T9l [TS] [26], acyclic matrices ^21j|, graded 
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positive-definite matrices [2D], scaled diagonally dominant matrices [3], totally positive matrices 
|31j . certain indefinite matrices |36| . and Cauchy matrices (as well as, more generally, matrices with 
displacement rank one) [15J. For such matrices, recent algorithmic advances (see [24[ I25j ) make the 
cost of achieving high relative accuracy comparable to that of alternative (and less accurate) SVD 
methods. 

The con-eigenvalue algorithm considered here is based on computing the eigenvalue decomposition 
of the product, CC, of positive-definite Cauchy matrices C and C, and is similar to the algorithm in 
[17| for the generalized eigenvalue decomposition, as well as the algorithm in f23] for the product SVD 
decomposition. We also rely on the algorithm in [15] for computing, with high relative accuracy, 
the Cholesky decomposition (with complete pivoting) C — (PL) (PL)* of a positive-definite 
Cauchy matrix C. However, since we are interested in computing only con-eigenvalues of some 
approximate size S, we stop Demmel's Cholesky algorithm once the diagonal elements Da are small 
with respect to S and the desired precision. Since the diagonal elements Da decay exponentially fast, 
this allows us to accurately compute con-eigenvalues of size S (and the associated con-eigenvectors) 

m o(n (log 5-1)') operations. We also modify the Cholesky decomposition algorithm in [l^ to 

yield high relative accuracy for Cauchy matrices dj — ^/ch\/^/ (1 — lijj), with 7^ exp(— Tj), 
where the real parts of the exponents, TZe{Tj), may be extremely small in magnitude. We observe 
that the error bounds developed in [23] are not applicable to our problem since the condition number 
of a Cauchy matrix cannot be appreciably reduced by scaling the rows and columns. In contrast, the 
error bounds developed in this paper yield high relative accuracy for all the computed con-eigenvalues 
larger than 5 (and high accuracy for the con-eigenvectors), as long as the n leading principal minors 
of L^L are well-conditioned, and the relative gap between the con-eigenvalues is not too small 
(we have always observed this to hold in practice). In particular, if 5 is chosen small enough, the 
full con-eigenvalue decomposition is obtained with high relative accuracy. The derivation of our 
error bounds makes crucial use of the component-wise perturbation theory developed in [20J for the 
singular vectors of graded matrices (see also [Mj), as well as the component-wise error analysis in 
[20] and f33] for the one-sided Jacobi method. We also use the error analysis given in [28] for the 
Householder QR method. We note that although our error estimates are much more pessimistic 
than what we observe in practice, they provide a framework for understanding the high accuracy of 
the con-eigenvalue algorithm of this paper. 

It has been an established practice, in both numerical analysis and signal processing, to use 
i^-type methods for representing functions. On the other hand, it has been understood for some 
time that nonlinear approximations may be far superior in achieving high accuracy with a minimal 
number of terms (see e.g., t2^)- However, in spite of many interesting results (see e.g., [321 l38l 
IT4l [39l l40l mi m [HI [22]), the widespread use of nonlinear approximations has been limited by a 
lack of efficient and accurate algorithms for computing them (particularly for functions with sharp 
changes or singularities) . Our algorithms provide the necessary tools for computing optimal nonlinear 
approximations via rational functions, and come with guaranteed accuracy bounds. We believe 
that these new accurate algorithms may greatly extend the practical use of L°° approximations in 
numerical analysis (see [27,) and signal processing (see [5J). 

In Section [2.11 we describe the reduction problem for rational functions, and connect its solution 
to a con-eigenvalue problem for positive definite Cauchy matrices. We then present new algorithms 
for solving the con-eigenvalue problem with high relative accuracy. We follow up in Section [3] with 
examples of using the reduction algorithm to construct and use optimal rational approximations 
for functions with singularities and sharp transitions. In Section [4] we verify the accuracy of the 
con-eigenvalue algorithm by comparing the con-eigenvalue decomposition of randomly generated 
Cauchy matrices with that obtained via standard algorithms in extended precision. In Section [SJ we 
prove that the con-eigenvalue algorithm achieves high relative accuracy and that the con-eigenvalue 
decomposition is stable with respect to small perturbations of the parameters defining the Cauchy 
matrix. Finally, Section [6] compares the reduction algorithm of this paper with other algorithms in 
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the literature for constructing optimal rational approximations. For the convenience of the reader 
we also provide relevant background material in Section [T) The proof of a technical proposition may 
be found in Appendix. 

2. Accurate con-eigenvalue decomposition (an informal derivation) 

2.1. Constructing optimal rational approximations via a con-eigenvalue problem. In or- 
der to motivate our con-eigenvalue algorithm, let us explain how the accurate computation of small 
con-eigenvalues and associated con-eigenvectors allows us to construct optimal rational approxima- 
tions. 

We consider an algorithm to find a rational approximation r(e^'^*^) to /(e^'^*^) in (II. 2p with a 
specified number of poles and with a (nearly) optimally small error in the L°°-norm. The algorithm 
is based on a theorem of Adamyan, Arov, and Krein (referred to below as the AAK Theorem) [3J. 
We note that the formulation given below in terms of a con-eigenvalue problem is similar to the 
approach taken in [Tl] and [5]. 

Given a target accuracy S for the error in the L°°-norm, the steps for computing the rational 
approximaiit r{z), 

— ' z — rii ^ — ' 1 — riiZ 

1=1 i=l 

are as follows: 

(1) Compute a con-eigenvalue < Am < S and corresponding con-eigenvector u of the Cauchy 
matrix Cij = Cij(7j, a-j), 



(2.1) Cu = XmU, where u 



U2 



Cij 



V Un j 

and Oi — ^Joi/^i, bj — y^Wj, Xi — yj ~ ^t7- The con-eigenvalues of C are labeled in 
non-increasing order, Ai > A2 > • • • > A,i. 

(2) Find the (exactly) m zeros r\j in the unit disk of the function 

(2.2) H.) = fET^- 

Am 1 - liZ 

The fact that there are exactly m zeros in the unit disk, corresponding to the index m of 
the con-eigenvalue Am, is a consequence of the AAK theorem. The poles of r(z) are given 
by the zeros 77^ of v(z). 

(3) Find the residues /3m of r(z) by solving the to x to linear system 



The L°°-error of the resulting rational approximation r(e^'^") satisfies ||/ — t'IIq^ ~ Am, and is close 
to the best error in the L°°-norm achievable by rational functions with no more than to poles in 
the unit disk. Hence, we are led to the problem of computing, to high relative accuracy, small 
con-eigenvalues and the associated con-eigenvectors of positive-definite Cauchy matrices. 

In many applications it is natural (and advisable) to maintain the poles 7j in the form 7^ — 
exp(— Tj) (see e.g., [SI [8]). As we explain in Section[3l this is particularly important if the function 
/(e^'^") has singularities or sharp transitions. The advantage of this form is that, on a logarithmic 
scale, the nodes are well separated. In such cases, our algorithm computes the new poles 77^ — 
exp (— Ci) with nearly full precision in the exponents i.e., Ci ~ / ICd is close to machine precision 
even if is close to zero. 
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Remark 1. In practice, finding the new poles rji using the formula for v{z) in (12. 2p is ill-advised, 
since evaluating v(z) in this form could result in loss of significant digits through catastrophic 
cancellation. Indeed, it turns out (see (6, Section 6] and ^7\) that the values of the con-eigenvector 
components satisfy Ui — y/alv {--^i), i = It then follows that the sum (|2.2I) must suffer 

cancellation of about logj^o (-^m^) digits if v(7i) and v {z) are of comparable size (note that Am 
controls the approximation error and, thus, is necessarily small). On the other hand, the function 
values V (7^) — Uij ^/al, i — 1, . . . ,n, along with the n poles of v{z), completely determine (|2.2p . 
Since the poles 7^ of f{z) are often close to the poles rji of r{z), we have observed that evaluating 
v{z) by using rational interpolation via continued fractions with the known values v (7^) allows us 
to obtain the new poles r]i with nearly full precision. In particular, an approximation v{z) to v(z) is 
computed via continued fractions, 

(2.4) v{z) 



1 + 02 (z - 71) / (1 + 03 (z - 72) / (1 + •••)) ' 

where the coefficients Uj are determined from the interpolation conditions v{'yi) = v {'ji). If the poles 
7i are given in the form 7^ = exp(— t^), we find that Newton's method on w(exp(— ?y)) yields the 
new poles r]i — exp (— Ci) with nearly full relative accuracy even when Re (Ci) ^ 1; see Section |3] for 
more details (achieving high relative accuracy also requires slightly modifying the recursion formulas 
for the continued fraction coefficients a^). A more detailed description of the root-finding algorithm 
may be found in |27| . 

2.2. Accurate con-eigenvalue decompositions of positive-definite matrices with RRDs. 

The con-eigenvalue problem for a positive-definite Cauchy matrix Cij — aibj / {xi + yj) reduces to 
an eigenvalue problem, 

(2.5) CCu = \Cu = \\\^ u. 

We first discuss a somewhat more general problem of computing accurate eigenvalues and eigenvec- 
tors of matrices of the form AA^ where we assume that A has a factorization A = XD^X*, with 
X a (well-conditioned) n x m matrix (to < n) and 13 an to x m diagonal matrix with positive, 
non-increasing diagonal entries. The rectangular form of the factorization, m < n, will be important 
in the sequel. 

Let us define the toxto matrix D (X'^X) D, and consider its SVD, G = W^V*. Then G*G = 
VT,^V*, and the ith right singular vector (1 < i < to), = V{:, i), satisfies [DX*liD) (DX^XD) v, 
E^jW^. It then follows that Zi = XDvi is an eigenvector of AA with eigenvalue Ef^, since 

AAz^ = {XD'^X*) (XD'^X^) z, = 

= XD {DX*XD) (DX^XD) V, = ^iXDn, = Y.\zi. 



and, thus, Zi = XDvi is an eigenvector of AA. To summarize: given the decomposition A = XD X*, 
an eigenvector Zi (i < m) of AA is given hy zl = X (Duil^^j^^^j, where Vi is the ith right singular 



vector of the rnxm matrix G = D (^X^X) D. Here E^^ is the ith singular value of G, and the ith con- 
eigenvalue of A. Let us now present an algorithm for accurately computing the con-eigenvalues and 
con-eigenvectors of A (its derivation also relies on the background material collected in Section [7]) . 
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Algorithm 1 ConEig_RRD (X, _D) computes accurate con-eigenvalue decomposition of XD^X*. 
Input: rank-revealing factors X and D (of dimensions n x m and m x m), where the diagonal of 
_D > is decreasing. Output: m con-eigenvalues/con-eigenvectors of XDX*, contained in E and T. 

(E, T) ^ ConEig_RRD {X, D) 

1 . Form G = D {X'^X) D 

2. Compute QR factors (Q, i?) -S- Householder_QR of G (G = 
QK) , with optional pivoting (see Section 17 . 3D 

3. Compute the SVD factors ([/;, E, f/^) Jacobi (i?) of R (R = 

UiT,U*) , using one-sided Jacobi, applied from the left (see Section 17 . 4D 

4. Compute i?i = D-^RD-\ Xi = D-^UiY}/"^ , and Yi = 
R^^Xi (see ( [276D below) 

5. Form the matrix of con-eigenvectors T = 

XYi, and output con-eigenvalues E and con-eigenvectors T 



Importantly, for Cauchy matrices {A — C) the elements of D decay exponentially fast, and it 

1/2 

would appear that computing the con-eigenvectors "zl = XDUifEj might lead to wildly inaccurate 
results even if the right singular vector of G, Vi, is computed accurately. However, as we show 
in Section \E\ Algorithm [1] achieves high accuracy despite the extreme ill-conditioning of D. The 
key reason is that the right singular vector Vi, corresponding to the singular value E^i, scales like 

(1/2 1/2 \ ^ 

Djj/Ti^l , Ej/ /Djj 1 , and the computed singular vector Vi is accurate relative to 

the scaling in D and E in the sense that 

\v. (j) - V, < min ^1 O (6) . 

(1/2 1/2 \ 
Djj/Tij^l 1 S,;/ /Djj) decreases exponentially fast away 

from the diagonal i — j. 

Let us give an informal explanation of the reasons why Algorithm [1] yields accurate results. 
As discussed in Section 17.31 the QR Householder algorithm computes an accurate rank-revealing 
decomposition of G = QR. It turns out (see Lemma [T^ that R may be factored as i? = D^Rq, 
where Rq is graded relative to D in the sense that||Z?_Ro-D^^ || and ||Z?i?(7^I?~^ || are not too large, as 
long as the n leading principal minors of X^X are well-conditioned. Therefore, from the discussion 
in Section [7j3] (see in particular Theorem [TS)) . the one-sided Jacobi algorithm computes the ith left 

{1/2 1/2 1 
Djj/Yi^l , Ej/ /Djj>. It follows that 

D~^UiY]\['^ may also be computed accurately. Finally, since the ith right singular vector Vi of R 
(and G) satisfies 

(2.6) = {DR^D-'Y' (D-'uj:\i^) , 

the con-eigenvector H = X (^DvlYiY^^^ niay be computed accurately, as long as DRqD^^ is com- 
puted accurately and is well-conditioned (we show this is the case if n leading principal minors of 
X^ X are well-conditioned). The last step in Algorithm [T] uses the approach in [25J for computing 
highly accurate right singular vectors via solving a triangular linear system of equations. 

Remark 2. To obtain optimal rational approximations (see Section [2.ip . we need to compute small 
con-eigenvalues (and the associated con-eigenvectors) of Cauchy matrices of the slightly different 
form, Ci,j = ^/a'i^/&~j/ (1 - ^i^), i.e., with ai = ^Jalj^i, bj = = ^ ~T7- "^^'^ 

same reasoning as in [15] shows that the Cholesky computation of G (see Section 17. 2|) is performed 
with high relative accuracy, as long as the differences 7~ — % are computed with high relative 
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accuracy. As explained in the next section, 7^ — 7; may be accurately computed if 7^ is of the form 
7i — exp (— Tj), where the exponents are known accurately (see Section [3] for examples). 



Remark 3. Computing the normalized eigenvector u via (12. 5p determines the con-eigenvector, the 
solution of p.ip . only up to an unknown phase factor e"*"^/^. Indeed, given any solution A and u of 
(|2.5p and an arbitrary phase factor e"*"^, it is easy to see that Ae~*"^ and ite"*"^/^ also satisfy (|2.ip . 
Let us now determine the phase so that the con-eigenvalues A are positive. To do so, we compute 
the usual inner product (C (ue-**/^) , ue-*'^/^) = A (ue*'^/^, ue-*'*/^) and choose (p so that A > 0. 
Since C is a positive-definite matrix, it follows that (ue*"^/^, ue"*"^/^) > 0. From this we obtain the 
phase factor as e*"^ = / 



2.3. Accurate con-eigenvalue decompositions of positive-definite Cauchy matrices. If 

A = C is a positive-definite Cauchy matrix, then the modified GECP algorithm in ^15J computes the 
Cholesky decomposition C — (PL) {PL)* with high relative accuracy (see Section m|) . There- 
fore, Algorithm [T] for the eigenvalue problem of CC may be used, with X = PL, to compute all the 
eigenvalues and eigenvectors (and, therefore, the con-eigenvectors and con-eigenvalues of C). 

For our purposes, we are only interested in computing a single con-eigenvector with associated 
con-eigenvalue of approximate size S (see Section [O]) . However, the diagonal elements of D may 
be many orders of magnitude smaller than S, and it is then natural to expect that, by computing 
a partial Cholesky decomposition of C, we may obtain the ith con-eigenvector in much fewer than 
O (n'^) operations. In this case, we stop Demmel's algorithm for the Cholesky decomposition of C 
once the diagonal elements Df^ are small with respect to the product of and the machine round- 
off e, that is, as soon as < 5^e for some m (notice that complete pivoting ensures that the 

diagonal elements Da are non-increasing). We then obtain C ~ C — (P^j (^^) ' ^^^^^ ^ 

an TO X n matrix, L is an n x m matrix and _D is a diagonal m x m matrix. Algorithms [2] and [3] 
contain pseudo-code for computing L, D, and P. In the pseudo-code / (n, to) denotes the first m < n 
columns of the n x n identity matrix. 
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Algorithm 2 Pivot _Order {a, b,x,y, 6) pre-computes pivot order for Cholesky factorization oi n x 
n positive-definite Cauchy matrix Cij = aibj/ (xi + yj). Input: a, b, x, and y defining = 
aibj/ (xi + and target size S of con-eigenvalue. Output: correctly pivoted vectors a, b, x, and y, 
truncation size m, and m x n permutation matrix P 

^a, 6, X, y, P, rnj -s— Pivot_Order (a, b, x, y, 6) 

Form vector gi := aibi/ {xi -\- yi) , i = 1, . . . ,n 
Set cutoff for GECP termination: 77 := eS^ 

Initialize permutation matrix (n x n identity): P = I{n,n) 
Compute correctly pivoted vectors: 
m := 1 

while \g (m)| > rj or m = n — I 

Find m < I < n such that \g{l)\ = max|g(m : n)\ 
Swap elements: 
g{l)^g{m), x(l) ^ x{m) , y{l) ^ y{m) 
a{l) O a{m) ,b{l) O b{m) 
Swap rows of permutation matrix: 

P(/,:)oP(m,:) 
Update diagonal of Schur complement: 
g{m -\- 1 : n) :— 
{x (m + 1 : n) — x{m)) / {y (m -|- 1 : n) — y{m)) g{m -f- 1 : n) 
Increment iteration count: 
m := m + 1 
Output a, 6, x, y,P {\ : rriju) ,m 



Algorithm 3 Cholesky_Cauchy (z, y, a, 6, (5) computes partial Cholesky factorization of positive- 
definite Cauchy matrix dj = aibj/ {xi +yj)- Input: a, b, x, and y defining dj = aibj/ {xi +yj), 
and target size S of con-eigenvalue. Output: n x m matrix L, m x m matrix Z), and permutation 
m X n matrix P in partial Cholesky factorization. 

(^L, D, P^ -i— Cholesky_Cauchy (a, 6, x, y, S) 
Compute pivoted vectors and matrix size m (Algorithm [5P : 

^a, b, X, y, P, rnj Pivot_Order (a, 6, x, y, S) 
Initialize generators: 

a := a, P :=b 
Compute first column of Schur complement: 
G(:,l) :=«*/?/ {x + y) 
for k = 2,11X1 

Update generators: 

a{k : n) :— a{k : n) * {x {k : n) — x {k — 1)) / (x {k : n) + y (k — 1)) 
I3{k:n):^l3{k:n)*{y{k:n)-y{k~ 1)) / (y (fc : n) x (fc - 1)) 
Extract fcth column for Cholesky factors: 

G {k : n, k) := a{k : n) * /3 {k : n) / (x {k : n) + y (k : n)) 
Output partial Cholesky factors: 

D = diag (G(l : n,\ : mfl'^ , Z = tril (G(l : 1 : to)) I)-^ + / (n, m) , P 
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Once the partial Cholesky decomposition C ~ C = yP^j \J computed, Algorithm[l]for 

the eigenvalue problem of CC may then be used, with X = PL and D = D, to compute accurate con- 
eigenvalues and con-eigenvectors of C (see Theorem[71) . Since the con-eigenvalues decay exponentially 
fast, the complexity of this algorithm is O (n(Jog{5e)~^)^^ operations. Therefore, when used in 
the reduction procedure outlined in Section 12.11 the near optimal rational approximation may be 
obtained by computing the SVD of a matrix that is roughly twice the size of the optimal number of 
poles. The pseudo-code is given in Algorithm |4l 



Algorithm 4 Con_Eigvector (a, b, x, y, 6) computes accurate con-eigenvalue decomposition of 
positive-definite Cauchy matrix Cij = aibj / (xi + yj) . Input: a, 6, x, and y defining = 
aibj / {xi + yj), and target size S of con-eigenvalue. Output: con-eigenvalues lager than 6, and 
associated con-eigenvectors. 

(S, T) ^ Con_Eigvector (a, b, x, y, S) 

1. Compute partial Cholesky factors {L,D,P)^ 
Cholesky_Cauchy (a, 6, X, y, (5) (Algorithm [3]) and set X ~ PL 

2. Compute con-eigenvalues and con-eigenvectors (E,r)-s— 
ConEig_RRD ( A, £)) using Algorithm [1] 

3. Select largest I such that En > 

S and output T, {1 : 1,1 : I) , T (1 : n, 1 : /) 



Remark 4. In applications involving functions / (e^'^*^) with singularities or sharp transitions, the 
poles 7i are given in the form 7^ = exp (— t^), where TZerj > and < IrriTj < 27r and the exponents 
Ti are known with high relative accuracy. Indeed, this form naturally arises either via a discretization 
of an integral (see [6l [8]) or as a result of an intermediate computation as in [27]. This leads us 
to modify Algorithms [5] and [3] so that high relative accuracy is achieved for poles of this form. In 
particular, we modify formulas (|7.7p . (|7.8p and ()7.9p in Section [T] For example, the formula for af^ 
in (|7.9p involves computing 

Xj - Xk~i _ 7j ^ - 7fc-i _ l-exp(-Tj +Tk-i) 
Xj + yk-i ^J^ - tFT 1 - cxp (-Tj - TfeTT) ' 

The simple modification is to use the Taylor expansion 1 — exp (z) z + z^/2 4- . . . if is small. 
The other formulas in (|7.7p . (17. 8p and (|7.9p are modified in a similar fashion, allowing the LDU 
factorization of C to be computed with high relative accuracy. 

In Section [3J we consider a case where the absolute values of many poles agree with 1 to twelve 
digits (i.e., the poles 7^ satisfy |7j| « 0.999999999999raxx). 

3. Examples of optimal rational approximations 
In this section, we consider some applications of the reduction algorithm. 

3.1. Optimal rational approximations of functions with singularities. Using the reduction 
algorithm, as well as tools developed in [S', 'S], we construct a (near) optimal rational approximation 
of a (piecewise smooth) function / with a finite number of isolated integrable singularities. For 
simplicity, we assume that singularities of / are at two points, and xq. 

Performing integration by parts L times on the expression for the Fourier coefficients, 

/„ = f /(:r)e2— = H f{x)e^"^^dx + f f{x)e^^'"^dx, 

Jo Jo Jxo 

we obtain 
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(27rm) Jo 



2'Kinx 



dx 



(27rm) 



27rinx 



where 



'^'^ = E 7^-^ (e2™-«i^(f-i) (xo) + F(P-i) (0)) 
^ izmn) V / 



^ (2Triny 

a;^, X indicate directional limits. As the first step in construct- 
ing a (near) optimal rational approximation to /, we subtract the leading L terms of the asymptotic 
expansion of /„ and consider (7„ = fn — hn- Since g„ decays like O (l/n^^^), it is sufficient to use 
the algorithm in [6l [8] to construct an approximation 

M 



(3.1) 



E 

m— 1 



< e, n>l. 



This algorithm requires quadruple precision for computing small singular values of a Hankel matrix 
but, due to the fast decay of gn, the matrix is small so that the computational cost is insignificant. 
An alternative method for obtaining p.ip based on rational representations of B-splines requires 
only double precision and will appear elsewhere . For hn we use a discretization of the integral 
representation for l/nP in [5] to obtain 

A/2 



(3.2) 
where r, 



1 



< e, l<p<L, 1 < n, 



' '^™.P ~ (p-i) 
imply that there are at most O ^(loge^^ 

e, for all n > 1. Note that when m < the nodes 7^ 
Thus, we arrive at 

A/2 

hn- «me-(^"+2--")" 
m=-A/i 



'^Yji-e^'"" and h is the step size used in the discretization. Results in [8] 
terms in the approximation of for a given accuracy 
e~'^ w 1 — e''™ are very close to one. 



(3.3) 
where 



A/2 

E 

m=-A/i 



< 2e, 



1 _ 1 _ 

^ w o ^-^^^ (2^0) Im.p, — / ^ ^-^^^ ^ (0) am,p- 



p=l 

Combining the approximations p.ip and (|3.3p . we obtain the suboptimal approximation 

M M2 A/2 



(3.4) 



< 3e, 



m— 1 m— — A/i m— — A/i 

where the number of terms is excessive (for the accuracy 3e) . We now use the reduction algorithm on 
p.4p to obtain a nearly optimal number of terms to approximate the Fourier coefScients /„ for n > 1. 
This, in turn leads to a near optimal rational approximation to f{x) with a nearly equioscillating 
error. 

As an example, we apply this procedure to 



(3.5) 



sin(4/37ra;), < a; < 3/4 
^0 3/4<x<l 

Choosing the parameters Mi — 200, M2 — 10, and h — .316707 in p.4p (see [8] for how to select the 
parameters) yields a sub-optimal approximation containing 426 pairs of conjugate-reciprocal poles 
7j = e~'^^ , which approximates / {x) in the L°° norm with error « 5 x 10"^^. We note that many 
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-14.5 - 



0.2 0.4 0.6 0.8 1.0 - I 

(a) (b) 

Figure 3.1. (a) Error of the rational approximation to / (x) in (|3.5p . (b) A zoom 
on a neighbourhood around one of the sing ularities x (E (3/4- 10"l^ 3/4 + lO^^^) . 

of the poles are extremely close to the unit disk (the magnitudes |7i| sa .999999999999a;a;xx of over 
a dozen poles agree with 1 to twelve digits). 

We apply the reduction algorithm using the approximation error S = IQ-" (thus, the Cholesky 
decomposition algorithm [3] is truncated once the diagonal elements are smaller than eS^, where e 
denotes the machine roundoff). As explained in Remark^ Algorithms [5] and |3] are modified to accu- 
rately compute the partial Cholesky decomposition for poles in the form = e~'^^ . After applying 
the reduction algorithm with approximation error 5 = 10"^'^, the resulting rational approximation 
contains 92 pairs of conjugate-reciprocal poles (i.e., about 46 poles per singularity). The resulting 
error is shown in Figure 13.11 

We note that the only step of the reduction procedure where quadruple precision is used is in 
computing the residues /3j (see Step 3 of Section However, using the techniques described in the 
background Section [7^ to factor the m x m Cauchy matrix, this step takes only O (m^) operations, 
and so does not impact the overall speed of the algorithm (recall that m denotes the number of 
reduced poles). 

We find that the exponents, 7]i, of the near optimal poles C,i = exp {—rji) are computed with high 
relative accuracy, i.e., 

|Re(?7,) - Re(r7j)| < |Re (?7,)| \r], - t],\ < \i],\S2, 

where 6i < 1.48 x 10"^^ and 62 < 14.87 x 10"^^. As a gauge we used the poles Q obtained in 
Mathematica -^^^ via extended precision arithmetic. We note that the real parts of some of the 
exponents 77, are of size |Re (77^)! w 10^^^. 

3.2. Solving viscous Burgers' equation. In ^17] we use the reduction algorithm to solve viscous 
Burgers' equation, 

(3.6) Ut — uux = vuxx, u{x,0) — uq{x), u{0,t) = u{l,t), .Te[0, 1], t>0. 

The solution of this equation develops a shock (or a sharp transition) on an interval of size O {i'). 
We approximate solutions to (|3.6p using rational functions of the form 

-A-fo /,\ Jl-fo 777 

fr[ - 7j it) j-[ e2-- - 7, {t) 

The key idea is to develop a numerical calculus using the reduction algorithm. Although opera- 
tors such as multiplication and convolution increase the number of poles in the representation, the 
reduction algorithm is employed at each stage to keep the number of poles near optimally small. 
Overall, about 10^ applications of the reduction algorithm were employed to compute the solutions 
illustrated below, thus confirming its robustness and efficiency. 
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Figure shows the computed solutions u{x,htj) to (|3.6p . with the viscosity ly = 10~^, the 
step size ht and the initial condition uo{x) = sin(27ra;) + l/2sin(47ra;). In our reduction procedure, 
we used the step size of ht — 10^^ and the error tolerance S — 10^^ (to match the error of our 
time discretization). The solution u{x,htj) is shown for time steps tj = hfj, j — 10^, 10'', 2 x 
10'', 3 X lO'*, 5 X lO''. We see that the solution u{x,t) develops two moving sharp transition regions, 
which approach each other and eventually merge into a single one about x ~ 1/2. The rational 
representations of u{x,tj) have 4, 11, 33, 29, and 19 pairs of conjugate-reciprocal poles, respectively. 
It also demonstrates that the transition regions of u{x,t) occur within intervals of width of O (v). 




Figure 3.2. (a) Solution it(x,tj), for 10"^, .1, .2, .3, and .5. (b) u{x,tj) in the 
transition region (1/2 - 10-^ 1/2 + 10"^), for = 0.4 (from [27\). These solutions 
are represented with 4, 11, 33, 29, and 19 pairs of conjugate- reciprocal poles. 



4. Accuracy verification 



We test the accuracy of Algorithm |3] on 500 random Cauchy matrices, Cij = {ctiOj) / (1 — TiTi)' 
i, j = 1, . . . , 120. The complex poles jj — pje'^"'^^ and residues aj — (jc'^"'^^ are generated by taking 
Pj, 4>j, and ^jj from the uniform distribution on (0, 1), and taking from the uniform distribution 
on (0, 10). For each randomly generated matrix, we first compute, as a gauge, CC = Z11Z~^ using 
the in-built Mathematica eigenvalue solver with 300 digits of precision, and compare the result 
with Z and E computed via Algorithm 2] using standard double precision. We then evaluate the 
maximum relative error in the con-eigenvalues \j — Tijj, 

error in the computed con-eigenvectors, maxj 

by the complex-valued constant Z (io,j) jZ {iQ,j), io — niaxi<i<„ \Z (i,j)\, since Z and Z 
are defined only up to an arbitrary complex- valued factor. 

Figures |4T] and l4?2l summarize the result of a typical run. Figure l4TT a) shows the distribution of 
the poles inside the unit disk and Figure HTT b) displays logj^Q A| as a function of the index j. Fig- 

ures l4.2r a) I4.2r b) show the relative errors in the con-eigenvalues A, — A, 



Z{:,j)-Zi:,j) J\\Z{:,3 



Xj — A j / I Xj I , and the maximum 
2. We first scale Z (:, j) 



con-eigenvectors \\zj 



•^J 1 1 2 . 



/ ll^jlL, both as functions of the index j. 



I \ Xj \ and the normalized 



In Figures 14731 and l4~4l for each of the 500 random Cauchy matrices, we plot the error in the com- 
puted con-eigenvalues Xj — Xj I \Xj\ and con-eigenvectors ||£, — Zj||2 / HzjUjfor j ~ 1,40,80,120 

(note the exponential decay of A^). We see that the con-eigenvalues and the con-eigenvectors 
are computed with nearly full precision for all the Cauchy matrices. In fact, the largest errors 

At — At 



l\Xj\ and 



Zj|L / \ zj\„ in the computed con-eigenvalues and con-eigenvectors, for any 



of the 500 Cauchy matrices and any 1 < j < n, are 5.13 x 10 and 5.35 x 10 
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Plol of poles (jui's) 



(a) 




con-eigenvalues, log scale 




Figure 4.1. (a) Distribution of poles 7^ determining Cauchy matrix C in a typical 
run. (b) Exponential decay of the eigenvalues Xj of CC as a function of the index 
j using logj^Q scale. 

5. Accuracy and perturbation theory 

We show that Algorithm 0] of the previous section achieves high relative accuracy. We also 
demonstrate that small perturbations of Oi, bj, Xi, and yj determining C lead to small relative per- 
turbations of the con-eigenvalues and small perturbations of the angles between subspaces spanned 
by the con-eigenvectors, as long as the parameters Xi and yj are not too close in a relative sense. 

For two (complex) floating point numbers x and y, we denote by fl (a; y) the result of applying 
the operation x Q y in floating point arithmetic, where is one of the four basic operations, € 
— , X , -;-}. We use that & {xQy) = {xQ y) (1 + S), where |^| < ce + O (e^) , e denote the machine 
round-off, and c is a small constant (cf. [5^). We will also abuse notation by letting fl (XY) denote 
the result of multiplying matrices X and Y in floating point arithmetic. Finally, throughout this 
paper the notation ||-|| always denotes the Frobenius norm. 

In Theorems [5][7] below we always assume that the con-eigenvalues are simple, although this is not 
a crucial restriction. In the statements and proofs of these theorems, the implicit constant factor 
implied by the notation O (77) and O (e) (here e,ri ^ 1) depends only on the size n of the matrix C. 
We also use the notation O (1) to denote a quantity that depends only on the size n. We note that 
all these implicit constants may be tracked more carefully and are modest-sized functions of n. 

The bounds in the theorems below depend on the Cholesky factors in the decomposition C = 
(PL) (PL) . In particular, let B = L, and consider the (symmetric) LU factorization B = 
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(a) 



(b) 



error in con-eigenvalues, log scale 



20 • 40 60 



100 120 



error in con-eigenvectors, log scale 



Figure 4.2. (a) Relative error in the jth con-eigenvalue, Xj — Xj 

tion of the index j. (b) The error in the jth con-eigenvector, | 
Zj = Z as a function of the index j. 



/ I A j I , as a func- 



■J 1 1 2 ' 



LgL"^, where Lb is unit lower triangular. Then the estimates in Theorems [S] - [7] depend on the 
quantities 



(5.1) 



\Lj^^\\^ k{Lb) : 



^J■o (Lb) = 

Ml (is) = max|||Lj3^||^,||i_B||^|K(iB), 

M2 (is) = II^B^lf A*i (is) (is) , 

where the condition number k[Lb) = 1 1 is || ||ij3^|| is not too large if B and its n leading principal 
submatrices are well-conditioned. The estimates in Theorems |6][7] also depend on 

(5.2) fi3 {Lb) ^ ||i~i (pmV'M2 (is) + J^K^ (is)) , 

where p, n, and ip are "pivot growth" factors associated with the QR factorization (see Section [73)) . 
and the factor u is associated with the one-sided Jacobi algorithm (see (|7.13p ). 

Remark. There are simple formulas for Lij and (L^^)^^. f|10|) in terms of the parameters a^, bj, Xi 
and Hj defining the Cauchy matrix C , and it is possible that the bounds below may be improved by 
using this additional structure. 

Theorem 5. Suppose that the parameters defining the positive-definite Cauchy matrix C — C{a, 6, x, y) 
are perturbed to d — a + 5a, b = b + Sb, x = x + Sx, and y = y + 6y. Let us define 

ij = (l/m + l/V2 + l/m) niax{||^a|L , ll^6|L > , ll^ylLI : 

where 

. . \yt-yj\ . \x^ + yJ\ 
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(a) 



(b) 



(c) 



• V'»'i*v"T'..'»V'£;.\~,»-» .V." i-. I * • « . 



(d) 



300 400 500 

/|A,|, for 



Figure 4.3. Relative error in the computed con-eigenvalues, 
j = 1,40,80,120 ((a), (b), (c), and (d), respectively), plotted for each of the 500 
random Cauchy matrices. 

Let C = LDL* denote the Cholesky factorization ofC, and let C = C{a,b,x,y) denote the Cauchy 
matrix corresponding to the perturbed parameters. Finally, let Zi, z^ denote the con- eigenvectors of 
C and C , corresponding to con- eigenvalues \i and \i . 

Then the relative difference in the con-eigenvalues Xi and Xi is bounded as 



Xi — Xi 



X,, 



< no {LB)0{ri). 



and the acute angle between the con- eigenvectors zi and Zi is bounded by 

sin Zi) < K (L) f ^^^^ + MO (Lb) Mi {Lb)) O (t?) 



relgapi 

Here /io {Lb), Mi {Lb) o,nd M2 {Lb) are defined in Ii5.1\) . and 

relgap, = mm 

I Ail + \Xj\ 

Next we state 
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(c) 



(b) 



(d) 



Figure 4.4. Relative error in the computed con-eigenvectors, \\zj — Zj\\2/\\zj\\2, for 
j — 1,40,80,120 ((a), (b), (c), and (d), respectively), plotted for each of the 500 
random Cauchy matrices. 



Theorem 6. Suppose that Algorithm [7] is used to compute the full con- eigenvalue decomposition 
of a positive-definite Cauchy matrix C . Suppose also that C has the Cholesky factorization C — 
[PL)D^ (PL)* , where P is the permutation matrix that encodes complete pivoting. 

Then the relative error between the computed con- eigenvalue and the exact is bounded as 



A,; — A,; 



I A,: 



where p, p, and are "pivot growth" factors associated with the QR factorization (see Section ] 7. 3\} . 
and the factor v is associated with the one-sided Jacobi algorithm (see \7.13\) ). 

Letting Zi, Zi denote exact and computed con- eigenvectors of C , the acute angle between Zi and 
Zi then satisfies 

' P3 (Lb) 



sin {Zzi, Zi) < K (L) 



relgapi 



|i-if «3 {LB)]0{e), 



where relgap^ is defined as in Theorem\^ and /is {Lb) is defined in \5.S\ 



Theorem 7. Suppose Algorithm [7] is used to compute m approximate con- eigenvalues and con- 
eigenvectors of a positive-definite Cauchy matrix C. Suppose also that C has the Cholesky factor- 
ization C — (PL) (PL)* , where P is the permutation matrix that encodes complete pivoting. 
Assuming that < A,;e for some 1 < i < m, the following error bound holds for the computed 

con-eigenvalue Xi, 



Xi — A, 



< {pptPflo (Lb) + ^+ \\C\\ pI (Lb)) O (e) , 
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and the acute angle between Zi and Zi is bounded by 

sin (Zz., z.) < n (L) {^^^^^^^^^^^f^ + H^-i' {Lb)) O (e) . 

In the above estimates, p, fi, and are "pivot growth" factors associated with the QR factorization 
(see Section \7.S\} . and the factor v is associated with the one-sided Jacobi algorithm (see 117. 13\) ). 

Remark 8. We note that the constants in the theorems above are significantly more pessimistic than 
we actually observe in numerical experiments. Indeed, while the bounds on the con-eigenvectors de- 
pend only on the well-conditioned matrices k (Lb) and k (L) (and, in particular, are independent of 
the exponentially decaying diagonal matrix D), they still scale like , where k = max {k (Lb) , k (L)}; 
the bounds on the con-eigenvalues are better — they scale like k^. However, in practice Algorithm 3] 
achieves nearly full precision for all the con-eigenvalues and con-eigenvectors. While it is likely that 
better estimates can be obtained, those presented here elucidate the basic mechanism behind the 
high accuracy that we observe in our experiments. 

5.1. Perturbation theorem. In this section, we prove Theorem^ We start by formulating several 
preliminary results. Lemma [S] describes how perturbations of the vectors a, &, x, and y defining 
the Cauchy matrix C — C{a,b,x,y) change the factors L and D in the Cholesky decomposition 
C = LDL* (see [15] for a proof). 

Lemma 9. Suppose the data defining the Cauchy matrix C — C{a, b, x, y) is perturbed to a = a + Sa, 
b — b + 6b, X = X + Sx, and y = y + Sy. Let us define 

r) = (l/?7i + l/r/2 + 1/V3) max{||Ja||^ , , , ||^y||^} , 

where 

= min ~ = min 1^' ~ ^ niin 1^' 

i^j \xj \ + \xi\' i^j \yj\ + \yi\' i^'j \xi \ + \y.j\ ' 

Then C = C{a, b, x, y) and C = C'(a, b, x, y) have Cholesky factorizations C — LDL* and C — 
LDL* , where L, L are unit lower triangular matrices, D, D are diagonal matrices with positive 
entries, and 



Lij Lij 



= |L,j|0(r?), Du-D,, = \Du\Oi7^) 



We now state the main result needed to prove Theorem [Sj Proposition [TO] considers how a 
perturbation B + SB in DBD affects the singular vectors. It turns out that, if all the principal 
minors of B are well-conditioned, then the errors in the perturbed singular vectors are graded, 

\v^{J) - - niin i\\SB\\) . 



I, V ^JJ ) 

The proof uses techniques developed in [4], [20], and [34]. 

Proposition 10. Suppose that G — DBD and G = D (B + SB) D, where B is complex symmetric 
(B^ B'^) and non-singular, and D is a diagonal matrix with positive, decreasing diagonal elements. 
Assume also that B has the LU decomposition B — LbL]^. 
Then the ith singular values "En and Tin of G and G satisfy 

<t.oiLB)Oi\\SB\\) 

where /xq (Lb) is defined in 15. il The ith (left or right) singular vectors Ui and ui of G and G, 
corresponding to (simple) singular values Tin o-nd Ha, may be chosen so that 

t1/2 ■ 

\ui{j)\ < Ml (is) min' 



3J 
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and 



- < ^2 (Lb) min ■ 



1/2 ■ 



0(\\SB\ 



where fii {Lb) and /i2 (Lb) are defined in \5.1l Finally, we have the following norm-wise bounds, 

Dui 



\B 



1/2 



(5.3) 



where /iq {Lb) is defined in \5.1 



< 



.1/2 



1-1/2 



< 



Proof. See the Appendix for proofs of the component-wise bounds on the left singular vectors. The 



bound on the relative difference 



is proven in |17j . 



1/2 

We now prove the norm-wise bound on Du,/T.l[ . To do so, note that DBD has the SVD 
DBD = IfSU^ , since B is complex symmetric (i.e., B has a Takagi factorization). Now suppose 

1 /2 

that DBDui = Y^uUi. Then the bound for Dui/T,^l follows from 

S,, = \{DBDm,u^)\ = \{BDu-,Dui)\ < \\B\\ \\Du^\\\ 
Similarly, note that D^^ B^^ D^^Ui — E^^u7. The bound for D^^UjYj^^l"^ then follows from 



\{D-^B-^D-h 



We now prove Theorem [SJ 



{B-^D-^u,,D-^u-i) \ < \\B- 



□ 



Proof. Recall that the matrix Z of con-eigenvectors satisfies Z — PL {^DVT, ^/^) , where C = 
{PL) {PL)* is the Cholesky factorization of C (with complete pivoting) and V is the matrix of 



right singular vectors oi G = D {L'^L) D. Let C 
tion of C. From Lemma[S] (see also |15j). 



PL ) ( PL ) denote the Cholesky factoriza- 



(5.4) 



L- L 



\L\\0{v), 



D, 



o{v). 



Defining G = D (^'^i) D, the above bounds yield G = L» {L^ L + E) D, where ||£'|| = O (r/) (since 
complete pivoting is used in the Cholesky factorization of C, \\L\\ ~ 0(1)). Since the ith con- 



eigenvalue of G is given by = Tiu, the estimate for the relative difference 
follows from Proposition [TUl 

To derive bounds for the singular vectors, we apply Proposition [TU] to G — D {iJ- L + E) D. In 
particular, there exist right unit singular vectors Vi of G and Vi of G such that 



(5.5) 
and 

(5.6) 
where 



\v,{3)-v,{])\<d,,{D,ll) 



\v^{J)\<^ll {LB)d,,{D,j:), 

Ai2 {Lb) 



relgapj 



o(||£;||)<d.,(AS)#^o(^) 



relgapj 



dij {D, S) = min ■ 



yl/2' D 



.1/2 ■ 
jj 



Here /ii {Lb) and ^2 {Lb) are defined in (|5.ip . Proposition [TOl also shows that 



(5.7) 



MO {LB)0{rj)<^ <1 



^^o{LB)0{7J). 
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Finally, defining Wi — Dvi/T^^l"^ and Wi — Dvi/T,n ^ , we have 



D 



n 



U/2 



1/2 



< 



< 



a/2 



{\v^{])-V^{3)\ + \v^{J)\^i^ {L b) O {7^)) 



m [Lb) 



+ fio (is) Ml (Lb) ]0{t]), 



relgapi 

where we used ()5.4p and (|5.7p in the first inequality, and (|5.5p - ()5.6p in the last one. 

Proposition [TOl also implies that 1/ \\wi\\ < Umax (L^L) = O (1) and 1/ < O (1) (since com- 
plete pivoting is used, ||L|| — O (1)). The proof now follows upon noting that the con-eigenvectors 
Zi and Zi satisfy Zi = Lwi, Zi — Lwl, and using ([57 

\Lwi — Lwi 



< 



< 



\Lw, 



\L\\ 



1^ 



\L\\ L- 



< k{L) 



M2 (Lb) 
relgapi 



/io (is) Ail (Lb) ]Oir]). 



□ 



5.2. Proof of Theorem [6] (high relative accuracy of Algorithm |4|). We now show that Algo- 
rithm |4] accurately computes the eigenvectors of CC (recall that C = {PL) (PL)*), as long as the 
n leading principal minors of L'^L are well-conditioned and the relative gap between the eigenvalues 
is not too small. 

Before proving Theorem |6l we first need several lemmas on graded matrices. 

Lemma 11. Let D be a positive definite diagonal matrix with decreasing diagonal elements, and let 
L and R denote nonsingular lower and upper triangular matrices. Then 



and 



\DLD-^\ < \\L\ 



\D-'^RD\\ < \\R\\ 



{DLD- 



[d-^rdV^ 



< L~ 



< R' 



Proof. Since the diagonal elements of D are decreasing and L is lower triangular, 



(DLD-^) 



< 



\L^J\ < \L.,,\ 



and 



{dld- 



< 



< 



L- 



Since the Frobenius norm is absolute, the first two inequalities in Lemma [TT] follow. The other two 
inequalities can be shown in a similar fashion. □ 



Lemma 12. Let D > denote a diagonal matrix with decreasing diagonal elements, and let B 
denote a non-singular, complex symmetric matrix with LU factorization B — LsLg. Suppose that 
DBD has the QR factorization QR — DBD. Then the upper triangular matrix Rq = D^^R satisfies 

\\DRoD-H\ < WLef , \\DRg^D-^\ < Ili^Mf ■ 
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Moreover, the ith left singular vector ui of R (corresponding to singular value "Su) satisfies 



(5.8) 



> \\B-^ 



-1/2 



(Lb) 



Proof First note DBD = [DLbD^^) (DL^D). Now suppose that DLbD^^ has the QR factoriza- 
tion QiRi = DLbD-\ Then Qi (RiDL^iD) = QR. Since RiDL^D is upper triangular, there is 
a diagonal matrix ft such that {flijl ^ I, Q = Qifl^^, and R = ilRiDL^D. Therefore, we have 

Ro = D-^R = VlD-^RiDLlD. 

It follows that DRqD-^ = fl {D-^RiD) L^, and we obtain 



\DRoD-H 



\n{D-'RiD) Ll\\ < ||i?i|| IlL^II < \\Li 



In the first inequality, we use Lemma and in the last one we use 
Lemma [TT] Similarly, we have 



\DLrD~^\\ and 



\dr;;^d-^ 



L 





<\\Ri'\\ 




< 











where we use Lemma I 111 in the first inequality, and ||i?j^i|| = HZ^i^^ZJ-i || and Lemma 1111 in the 
last one. 

In order to prove the bound for D^^UiT,ll^ , we first claim that ||Z3-1(3_D|| < k{Lb)- Indeed, 
since QiRi — DLbD^^ and Q = Qiil^^, we have 

||L»-igi:>|| = \\d-'^QiD\\ = \\Lb {D-^RyD)\\ 

< II^bII ||i?ri < ll^sll \\DLs^D-'\\ < \\Lb\\ IlLsi . 

In the above string of inequalities, we use of Lemma [TT] repeatedly 

Now, if i? has the SVD R = UY.V*, then D-^R-^D-^ has the SVD D-^R-^D-^ = FE"! (QU)* 
(recall that QR = DBD). Therefore, the left and right singular vectors Ui and Vi of R satisfy 
{D~^B~^D~^^ {Qui) — y^j^Vi. Since B is complex symmetric, we may also assume (without loss of 

generality) that Qui - 



The bound on 



D-^u^Y. 



1/2 



now follows from 



(B-i (D-^v^) ,D-^v,)\ < \\B-^\ \\D-\,\\ \\D-^v- 



= \\B-m\D-\£ =\\B-^ 



\D-^Qu, 



= \\B-^\\\\{D-^QD) D-\,\\ 
< \\B-'^\\k^ {LB)\\D-\if . 

In the last inequality, we used the bound ||Z3-i(3I?|| < k (Lb), as shown in the previous paragraph. 

□ 



Lemma 13. Let DBD, Q, R, and Lb be as in Lemma and define Ri = D^^RD"^ . Then 
k{Ri) < k?{Lb), and the ith right singular vector Vi of R (corresponding to singular value Y^u) 
satisfies 

= R^ (D-\,Yir' 



where 



\Ri\\<\\L, 



\R-A\<\\LBt^ 
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Proof. Define Rq ^ D-'^R and Ri = D'^RD-^. Then since Ri = DRaD~^, tlie bounds for 
and follow from Lemma [T^ Using (|2.6p . we also have that 



-1/2 



dr-\,j:^^ 



1/2 



(dr;;'d-') (d'^u^!' 



(DRoD-') 

R-' {p-^uj:l 



/2 



□ 



We now prove Theorem [51 



Proof. First consider Step 1 of Algorithm |4l From [T5], the computed Cholesky factors D and L of 
C satisfy 



(5.9) 



D, 



L- L 



<Oie) 



We now examine the error in applying Algorithm [T] to compute the con-eigenvectors and con- 
eigenvalues. 

In Step 1 of Algorithm[Tl the computed matrix G satisfies G = {l(^D (l^L^ ^ D {L^ L + Eq) D, 
where ||i?o|| < O (e) (recall that complete pivoting is used, so that ||L|| = O (1)). 

In Step 2 of Algorithm (TJ a computed upper triangular factor i? of G is obtained using the 
Householder QR algorithm. By Theorem [T7] in Section 17.31 there is an orthogonal matrix Q such 
that 



(5.10) 



QR^D(L^L + E2)D, 



where E2 = Eq + Ei, \\Ei\\ < piJ.ipO{e), and p, /i, and ip are "pivot growth factors" described 
in Section 17.31 Now suppose that L'^L + E2 has an LU factorization L + E2 = LbUb. By 
Lemma [T2l Ro = D~^R satisfies 

|2 



(5.11) 

where we used that \\Lb 



\DR^^D-^\\ < \\Lg^\ 



0{e)., 



Lf 



0[e). 



0{e) and 

Step 3 of Algorithm [1] involves computing an approximate SVD R w UiYiUr using the modified 
one-sided Jacobi algorithm, applied from the left. Note that, from (|5.10p . if R has the (exact) SVD 
R = Ui'EU*, then QUi and Ur are the matrices of left and right singular vectors of D {iJ^ L + E2) D. 

Moreover, if "En denotes the exact singular value of R (and D (^L'^L + E2) D), then Proposition [TOl 

ensures that 

^ </^o {LB)Oi\\E2\\) < pp^jpo {LB)O(e). 



Now let Yiii denote the computed singular value of R obtained via the one-sided Jacobi algorithm. 
Then from Theorem [T51 and the equality R = D^Rq, we also have that 



< i^oO (e) , 



Combining the previous two inequalities yields the bound on the computed con-eigenvalues (recall 
that the exact and computed con-eigenvalues satisfy — Y^u and Xi = T^u). 
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Now let Ui and Tin denote the ith computed left singular vector and singular value of R. Then by 
TheoremlTSland the equality R — D^Rq, there is an exact left singular vector Uj-^-* of R, corresponding 
to singular value Tn, such that 



(5.12) 



< 



relgap, 



where v is described in Theorem [T51 in Section [7j3] (see, in particular, (|7.13p ). 

Now define Ri — D~^RD~^ and Ri ~ &(^D^^RD^^^. We show that the computed vector 

y, w ^L^^ (B^Im^E^/^) , obtained from Step 4 of Algorithm^ is close to ^ = R];^ (^ZJ-^wf ^E^^ 
In particular, Step 4 involves computing an approximation jji to the triangular system Riyf^ = a 



where Xj — fi 



and 



i?i = fl D-'RD-' =Ri+SRi, 



\SRi\ 
\\Ri\\ 



<0{e). 



In the above inequality for ||(5i?i||, we used (|5.9I) . We will also need the following expression for 
which follows from Lemma [T^ 



(5.13) 



Ri = DRoD-\ K{Ri)<n'^{LB) + 0{€). 



Now, recall that y^^'^ is the exact solution of Riyf^ — Xi, Xi — D ^m^'eI/^. Since yi is computed 

from the system Riyf^ = Xi by backsubstitution, there is a matrix 5R2 such that ||(5-R2|| / ||^2|| = 
O (e) and (i?i + 5R2) yi — Xi + {xi — Xi). Therefore, 





yi 









(5.14) 



< KiRi) 



\{xi - Xi)\\ 



0{e) 



+ Oie) , 



where we used (|5.13p in the last inequality. To bound \\{xi — Xi)\\, we compute that 



f, (j) = fl ( DJ^u. ij) Sl/' ) - DJ^Hi ij) (l + O (e)) 



U/2 



1/2 



D 



where (15. 9p is used in the second equality to bound Djj / Djj and Proposition [TU] is used in the last 
equality to bound Tn/Yin. We also have from Lemma \T2\ (with B = L"^L + E2, II-B2II ^ pi^ipO (e)) 
that 



\\x^\\-' = 
Therefore, by (|5.12p . we have 









Ikili 


n_i (l)?^l/2 





< 



v'^L II K (Lb) 
relgapj 



< \\L-H\K{LB) + 0{e). 



\L~s^\\ K{LB)0{e) 



'^^iLn) ]0{e). 



L- 
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It finally follows from (|5.14p and the above inequality that 

(v\\L-^\tc' {Lb) 



(5.15) 



relgapj 



+ \\L-^'rn^LB)\0{e). 



Now, from Lemma [T51 there exists a right singular vector v^^^ of R such that y^^^ = Dv';pYn^^'^ . 
Moreover, from (|5.10p . is also a right singular vector of D {L'^L + E2) D, II-E2II < p^tpO (e). 
Therefore, Proposition [TU] ensures that there is a right singular vector vi of G — D {L^ L) D such 
that 



(5.16) 



< ^^0(||i?2||)<^^#^0(.), 



relgapi 



We also have from Proposition [TO] that 
D^^WjEV^, it follows from (IETO)l that 



(5.17) 
(5.18) 



n-l (l)?^l/2 



relgap, 

< Therefore, letting 



(1) 
yi - yi 










y?" 








y?" 





where we use that 
we finally obtain 



,(1) 



relgap, 

jy,|| + O (e) and sT/^' = + O (e). Combining (jCTf)) and (jETO)) . 



L^^ll «'^(Lb) 0(e), 



where 



\\yi\\ V relgapj 

A*3 (-^s) = ll-^^^ll (^^"0^2 (is) + i^K''' (-^b)) ■ 

In the final step of the algorithm, we compute an approximation = fi (^Ljjij 

eigenvector Zi = Lyi. From (|5.9p . we have = LyJ + Cj, where ||ei|| < \\yi\\ O (e). Therefore, we 
obtain 



to the true con- 



< 



\Lyi - Lyi 

\\Lyz\\ 



\L\ 



< \\L\\ L- 



\\yz 



0{e) 



< k{L) 



M3 (Lb) 
relgapj 



iL^'fn' {LB)]0{e). 



□ 



5.3. Proof of Theorem [3 We only prove the error bounds for the computed con-eigenvector 
components in Theorem [7] (the error bounds for the computed con-eigenvalues follow in a similar 
fashion) . 

We need the following well-known result describing the sensitivity of the eigenvalue problem for 
diagonalizable matrices. 

Lemma 14. Assuming that A has simple eigenvalues, we consider its perturbation A E. Let 
Zq — [ zi ... Zn ) denote a matrix of unit eigenvectors of A, with corresponding eigenvalues 
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Ai,...,A„. Then, for small enough \\E\\, the ith unit eigenvectors Zi and Zi of A and A may he 
chosen so that 

Q (II -Ell ) 

-, where ahsgap^ — min |Ai — Aj| . 



\zi - ZiW < k{Zo) 



ahsgapj^ 



The next result shows that the matrix of eigenvectors of CC is well-conditioned. Recall that O (1) 
denotes a constant that depends on n only (C has dimensions n x n). 

Lemma 15. Let C denote a positive-definite Cauchy matrix, and let Zq denote the matrix of unit 
eigenvectors of CC . Then we have k{Zo) < k {L) {Lb) O {!) . 

Proof. From Section 12.21 we know that the matrix Z of (unnormalized) eigenvectors of CC is given 
by Z = PL (DVE^^/^), where V is the matrix of right singular vectors of D {L^ L) D and L is the 
lower triangular matrix in the Cholesky factorization C = (PL) {PL)* . Now define the matrix, 
Zq, of normalized eigenvectors Zq = Zn~^, where flu — \\zi\\. 
By Proposition [TOl we have that 



, < ^1 {Lb) 



where jii (Lb) is defined in (|5.1I) and denotes the matrix of eigenvalues of CC. Therefore, 
ll^^ll < y^ll^ll < V^^\\L\\^ll {Lb). Also, 



Z' 



:z-') 



^1/2 



V D- 



{PL)-' 



< 



L- 



= L- 



1/2 



< L- 



where we used that V =Vm. the last equality. It follows that 



Zn 



< lirill \\Z' 



< y^K {L)til {Lb). 



Finally, using the above inequality and the bound ||2'o|| < ^/n (recall that the column norms of Zq 
are unity), 

K {Zq) < UK {L) {Lb) ■ 

□ 

The next lemma is the key to proving Theorem [71 

Proposition 16. Suppose that Algorithm [5| produces, in exact arithemtic, the partial Cholesky 
factorization C = ( PL J I PL ) , where P has dimension m x n, L has dimension n x m, and 



D has dimension m x m. Also assume that < eSfj for some (simple) eigenvalue Sfj of CC 

(1 < i < m). 

Then the ith unit eigenvectors Zi and zi of CC and CC may be chosen so that 
11^.^ _ ^^11 ^ \\C\n{Lp^ALBl^ ^^^^^ ^ ^ 

relgapi j^i E^, 



-'n 



where /^i {Lb) is defined in i5.1\} . 

Proof. After m steps of Gaussian elimination with complete pivoting. 



P^CP = 









^2*1 ^ 




\ F21 


F22 ) 




TT"* 

^22 J 





^11^11 ^11^21 
F21F11 C22 



where G22 = F2iF^^ + F22F22 and IF22 (i, j)| < Fn (m, m), i, j = m + 1, . . . ,n (recall that C > 0). 
Now, lb ' ^ ^ nT 



-F^ii ^21 ) and 
P^CP = 



i^ii 
-F21 



( ^1*1 ^^2*1 



^ll-f^ll ^11^21 
-^"21^11 ^21^21 
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Therefore, since D„ 



and 



Fii {m,m), 



C-C 



< Dl„P (1) , 



cc~cc 



<\\C\\Dl„P{l) 



It follows from Lemmas [T3] and [T5] that the ith eigenvectors Zi and Zi of CC and CC may be chosen 
so that 

^ rK{Z)\\C\\Dl^. 



\zi — zA\ < max 



< 



L>2 



n(Z) 



\C\\ 



relgap, " " Ef 
In the last inequality, we used that (see |34j) 



< 



\C\\K{L)^ll (Lb) 
relgap, 



0{e). 



a ^ jj 
max — iri 



< max ■ 



l^ii — Sj-j- 1 



□ 



Finally, we are ready to prove Theorem [T] 

Proof. The proof of Theorem [S] shows that the ith computed unit eigenvector zl of CC, C = 
, is close to an exact unit eigenvector Zj, i.e., 

V relgapj y 
where ^^^{Lb) is defined in (j5.2l) . Also, Proposition 1161 implies that there is an exact unit unit 



eigenvector Zi of CC such that 



z,; < O e) . 

relgap. 



The claim follows from the above inequalities. □ 

6. Comparison with related approaches for constructing optimal rational 

approximations 

Numerical approaches for finding near optimal rational approximations originate in theoretical 
results of Adamyan, Arov, and Krein [Tl|21[3]- In particular, given a periodic function / (e^'^™) £ 
L°°(0, 1), AAK theory yields an optimal "rational-like" approximation tm (e^'^*^), 

(6-1) rM [z) = -J- -T T- -— , 101 < 1, 

(z - Cij • ■ • (2 - Cm) 

constructed from the left and right singular vectors corresponding to the Afth singular value, ctm; 
of the infinite Hankel matrix Hij — f {i + j — l^) , i, j = 1, 2, . . .. The numerator of vm (z) in (j6.ip 
is analytic in the unit disk. The approximation error satisfies 

max|/(e2-") -rM(e2^")| =aM, 

where the number of poles C,j in (j6.ip equals the index M of the singular value om (index counting 
starts from zero). Moreover, the L°°-norm approximation error is optimal among all functions of 
the form ((6T|l . 

In order to use AAK theory to compute near optimal rational approximations, standard numer- 
ical approaches compute singular vectors of a truncated Hankel matrix. The poles of the rational 
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approximation are obtained as roots of a polynomial whose coefficients are the entries of the sin- 
gular vector. Such approaches have a long history of their own and, in particular, let us mention 
the pioneering papers [38, 3^1 140). A recent version (incorporating additional ideas) can be found in 

Instead of truncating the Hankel matrix, the approach of this paper is based on the observation 
that it is always possible (see e.g. [51 [51 [51 to construct a sub-optimal rational approximation, 
i.e., an approximation with excessive number of poles for a desired accuracy. This leads us to special- 
ize AAK theory to proper rational functions / (e^'^'^), and to formulate the reduction problem (see 
Section l2.ll and ^ Section 6]). Importantly, this results in a con-eigenvalue problem of finite size 
and with no additional approximations. Moreover, this formulation allows us to develop a numerical 
calculus based on rational functions (numerical operations such as addition and multiplication in- 
crease the number of poles; the reduction algorithm is applied to keep their number near optimally 
small, see [27]). Early approaches of this type can be found in [32l fT4l |4T] : however, these algorithms 
may require extended precision for high accuracy and also scale cubically in the number of original 
poles. 

Comparing our approach with that in e.g. [21|, we make two observations. First, to justify the 
truncation of an infinite Hankel matrix, the Fourier coefficients have to decay below the desired 
accuracy of approximation. Thus, for functions that have sharp transitions (as in the example of 
Section 3.2) or singularities (as in the example of Section 3.1), where the Fourier coefficients decay 
slowly, this would require computing singular values of very large matrices. In the examples of 
Sections 3.1 and 3.2, Hankel matrices of size « 10^ x lO'^ and ~ 10^ x 10^ would be needed in order 
to attain a comparable accuracy. This approach would also require finding roots of polynomials with 

10^ and « 10^ coefficients, respectively. 

Our second observation is that using Hankel matrices may require extended precision arithmetic 
if high accuracy is desired, as is the case in examples of Sections 3.1 and 3.2. Indeed, existing SVD 
algorithms do not accurately compute small singular values of Hankel matrices. Also, the roots of 
high degree polynomials (determined at the SVD step) may be sensitive to perturbations in their 
coefficients. However, when limited to approximating smooth functions, these "truncated Hankel" 
methods can yield surprisingly high accuracy since the errors in the poles may be compensated by 
the residues. As far as we are aware, truncated Hankel methods for constructing optimal rational 
approximations for functions with singularities generally do not achieve approximation errors better 
than sa 10""*. In contrast, in Section 3.1 we show that the reduction algorithm approximates 
piecewise smooth functions with errors close to machine precision. 

We also note that the results in [27] (illustrated in Section 3.2) demonstrate an effective numerical 
calculus based on the reduction algorithm, capable of computing highly accurate solutions to viscous 
Burgers' equation for viscosity as small as 10~^. These solutions exhibit moving transitions regions 
of width 10""", and computing them with high accuracy over long time intervals is a nontrivial 
task for any numerical method. The con-eigenvalue algorithm of this paper is critical to the high 
accuracy and efficiency of this numerical calculus. 



7. Background on algorithms for high relative accuracy 

Here we provide necessary background on computing highly accurate SVDs, as well some error 
bounds that are needed for the analysis of the con-eigenvalue algorithm. Although the results we 
need in [20l [33l [ITl [34l [151 HI] are only stated there for real-valued matrices, they carry over to 
complex-valued matrices with minor modifications and are formulated as such. 

7.1. Accurate SVDs of matrices with rank- revealing decompositions. According to the 
usual perturbation theory for the SVD (see e.g. [12J), perturbations 6 A of a matrix A change the ith 
singular value Ci by 6ai and corresponding unit eigenvector Ui by 6ui, where (assuming for simplicity 
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that ai is simple), 

(7.1) l^cr.l/cri < liMII , < 4^^, absgap^ =min|CT, -CTjI/cTi. 

absgapj i=ij 

Therefore, small perturbations in the elements of A may lead to large relative changes in the small 
singular values and the associated singular vectors. Moreover, since standard algorithms compute an 
SVD of some nearby matrix A + SA, where \\SA\\ / \\A\\ = O (e), the perturbation bound (|7.1|) shows 
that the computed small singular values and corresponding singular vectors will be inaccurate. 

In contrast, the authors in [T7] show that, for many structured matrices, the ith singular value 
(7i <C (Ti and the associated singular vector are robust with respect to small perturbations of the 
matrix that preserve its underlying structure. The sensitivity is instead governed by the ith relative 
gap 

, . \<^i-^j\ 

relgapj — mm — . 

ai + aj 

More precisely, let us consider the class of matrices for which a rank-revealing decomposition A = 
XDY* is available and may be computed accurately. Here X and Y are nxm well-conditioned ma- 
trices and Z? is an m X m diagonal matrix that contains any possible ill-conditioning of A. As is shown 
in [TT], a perturbation of A = XDY* that is of the form A + SA^ {X + SX) {D + SD) (Y + SY)* , 
where 

(72) M^-O(e) m-0(.) M^-O(e) 

^ ' 11X11 \\Y\\ lA.I " 

changes the ith singular value Ui and associated left (or right) singular vector Ui by amounts 5ai 
and 5ui bounded by 

(7.3) M<„,ax(.(X),.(r))0(e), ||<5..||<=M^ll^O(.), 

Oi relgap, 

where n{X) = ||X|j ||^^|| and X^ denotes the pseudo-inverse of A. One reason this class of matrices 
is so useful is that Gaussian elimination with complete pivoting (GECP) (or simple modifications) 
computes accurate rank-revealing decompositions of many types of structured matrices (see |17j and 
^15j). Moreover, small perturbations of such matrices that preserve their underlying structure lead 
to small perturbations in the rank-revealing factors and, therefore, small relative perturbations of 
the singular values. 

Given the decomposition A — XDY* , it is shown in [171 Algorithm 3.1] that an SVD of A may 
be computed with high relative accuracy, and with about the same cost as standard, less accurate 
SVD algorithms for dense matrices. The key to this algorithm is the one-sided Jacobi algorithm 
(briefly reviewed in Section [7.4p . which, with an appropriate stopping criterion, accurately computes 
the SVD of matrices of the form DB, where D is diagonal (and typically highly ill-conditioned) and 
B is well-conditioned (see [2D] and [33]). In particular, the algorithm in pTl Algorithm 3.1] yields 
computed singular values ct, and left (or right) singular vectors ui that satisfy 

(7.4) ^i^^ <niax(K(X),K(y))0(e), 

/_ -N ,, ^ t, niax (k (X) . K (y)) , , 

7.5 u,-u, < ^ \ ^ ^' 0{e), 

relgap, 

7.2. LDU factorization of Cauchy matrices. In this section we review how a modification of 
GECP computes accurate rank- revealing decompositions of Cauchy matrices [15J. 

We describe Demmel's algorithm (see Algorithms 3 and 4 in [15] and Algorithm 2.5 in [gj) for 
computing an accurate rank-revealing decomposition of a n x n positive-definite Cauchy matrix 
Cij = ttibj/ {xi + Hj) (note that Demmel refers to such matrices as quasi-Cauchy). The algorithm is 
based on a modification of Gaussian elimination for computing, in O (rt^) operations, the Cholesky 
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factorization C = {PL) D (PD)* of a positive-definite Cauchy matrix (more generally, the algorithm 
computes an LDU factorization for an arbitrary Cauchy matrix in O (n^) operations). Here P is a 
permutation matrix, L is a unit lower triangular matrix, and D is a, diagonal matrix with positive 
diagonal elements. It is shown in [15] that, remarkably, the components of the LDU factors L, U, 
and D are computed to high relative accuracy. 



(7.6) 



< \L^J\ 



< c„ \Uij\ 



< c„ I A 



where c„ is a modest-sized function of n. The basic reason the algorithm achieves high relative 
accuracy is that the only operations involved are multiplication and division of floating point numbers 
(additions and subtractions in the algorithm involve only Xi and yj, which are assumed to be exact). 

We now review the basic idea behind the algorithm in |15j . First, ignoring pivoting for a moment, 
we assume that, after k steps of Gaussian elimination, the Cauchy matrix is transformed to the 
matrix 

(k-\-l) (k) 

The elements of the Schur complement G22 may be computed from those of G22 by using the 
recursion 




(7.7) G« = ^1^) yi^)G\l'\ ^,J = k + l,. 



Introducing pivoting, we observe that the matrix G'-'^-' may be obtained by applying Gaussian elim- 
ination to a Cauchy matrix G'*-') = C^'^^ {a'^^\b^'^\x'^^\y'^^^) , where o''^^ b^'^\ x^^^ and y'^^^ are 
permutations of a, 6, x and y corresponding to the row and column pivoting of G. As long as the 
vectors a, 6, x and y are permuted according to the pivoting of G*-*^-*, the recursive formula (|7.7p still 
holds. 

It is observed in [TS] that if G is positive-definite (and, therefore, only diagonal pivoting is needed), 
then the pivot order may be determined in advance in O (n^) operations by computing diag (G^*^^) 
from formula (j7.7p . Once the correct pivot order is known, we do not need to compute the entire 
Schur complement G^'^^' to extract the components of L and ?7, but only its /cth row and fcth column. 
Indeed, we may use Algorithm 2.5 in [9], which uses the displacement structure of G, to compute an 
accurate Cholesky decomposition in O (n^) operations. To see how, note that it easily follows from 
(|7.7p that the Schur complement of a Cauchy matrix is a Cauchy matrix. 



^(fe)o(fe) 

(7.8) gW(z,j) = ^^, z,j=A + l,...,n, 

Xi + yj 

where the parameters a^^^ and /J^-'^'' satisfy the recursion 

C? n\ (fc) ~ (fc-1) M Vi - yk n{k-l) . , , , 

(7.9) a = ■ al PI = \ PI , i^k + l,...,n. 

Xi +yk Vi + Xk 

Since the fcth column L (:, fc) may be extracted from G'*^' (:, fc), we therefore only require O (n) oper- 
ations at each step of Gaussian elimination to compute L (:, fc). Updating a^"* and also requires 
only O (n) operations. In Section 12.31 (see Algorithms [2] and |3]), we present an O (^n (log (5^^)^^ 
algorithm to compute con-eigenvalues greater than a user specified cutoff S and, as a result, yield- 
ing a fast algorithm for obtaining nearly optimal rational approximations. Once an accurate LDU 

factorization G ~ ("^^) ^ (^"^) available, an accurate SVD of G may be obtained using the 
algorithm in [ITl Algorithm 3.1]. 
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7.3. Rank-revealing decompositions of graded matrices. We also review how a variant of the 
QR Householder algorithm with complete pivoting computes accurate rank-revealing decompositions 
of graded matrices [28J. 

It is shown in [5S] that the Householder QR algorithm with complete pivoting may be used to 
compute a rank-revealing decomposition of a graded matrix of the form A — DiBD2. Here Di and 
D2 are diagonal matrices that account for the ill-conditioning of A. Recall that the Householder 
QR algorithm uses repeated applications of orthogonal matrices to reduce A to an upper-triangular 
matrix R. On the first step, the parameter (3i and the vector vi of the Householder reflection matrix 
Q^^^ = I — Pivivl are chosen so that 



f an ^ 

021 



Q 



(1) 







V a„i / V / 

Consequently, the first application of Q^^-' to A results in a matrix of the form 

"11 "12 ■ • ■ "in \ 



V 

X (n 







,(1) 



,(1) 



and. after n — 1 such 



This process is repeated on the (n — 1) x (n — 1) lower block 

L J 2<i,j<n 

steps, ^("^1) = Q*^"^^^ . . .Q^^^A = R, where R is upper triangular. In the version considered in 
|28j . the rows of A are first pre-sorted so that so that \\A (1, :)||^ > • • • > ||^ {n, OIloo- "^^^^ algorithm 
then proceeds as above, except that at each step, fc, column pivoting is performed to ensure that 
II^W (fc : n,k)\\^ > ■■ ■ > \\A'-'''^ {k : n,n) ll^. Letting Pi denote the row permutation matrix that 
pre-sorts the rows of A, and letting P2 denote the column permutation matrix corresponding to the 
column pivoting, the QR Householder algorithm produces the QR factorization P1AP2 = QR. 

Following [28J, we consider the error analysis of the Householder algorithm (without pivoting) 
applied to P1AP2, where Pi and P2 are chosen so that no column or row exchanges are necessary (e.g. 
the matrix A is pre-pivoted) . Assume that the matrix P1AP2 may be factored as P1AP2 = D1BD2, 
where Di and D2 are diagonal matrices, and that the Householder algorithm, applied to the row- 
scaled matrix C — DiB, produces intermediate matrices C^^^ with columns c'-'\ Finally, define the 
quantities p, /i, and ^ by 



(7.10) 



max 



maxj^/j 




maxj 


Cij\ 



n = max max 

fe j>k 



Jk) 



{k : m) 



Sk) 



(k : m) 



ijj = max 



max," 



\Ckj\ 



l<i<n 
i<k<r, 



max," 



The above quantities measure the extent to which the Householder algorithm preserves the scaling 
in the intermediate matrices A^'^^ , and are almost always small (this is analogous to the pivot growth 
factor in Gaussian elimination with row pivoting). It is shown in [28J that 

Theorem 17. Suppose that A is pre-pivoted, and the Householder algorithm is used to compute the 
upper triangular matrix R of the QR decomposition. Then there is an orthogonal matrix Q such that 
QR ~ Di {B + SB) D2, where SB satisfies 

\\5B\\<pi,fi\\B\\0{e), 

and p, p, and ip are defined in 117. 1U\) . 

In |28j Theorem [17] is combined with the theory developed in fVl] (e.g., see Theorems 4.1 and 4.2 
in [17]) to show that the QR algorithm with complete pivoting produces accurate rank revealing 
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decompositions of graded matrices of the form A = D1BD2; as long as the principal minors of 
B are well-conditioned and the diagonal elements of Di and D2 are approximately decreasing in 
magnitude. 

Remark. Instead of pre-sorting the rows of A and applying the Householder algorithm with column 
pivoting, one may also use a version of the Householder algorithm in which both row and column 
pivoting is employed (see |28) for more details). Gaussian elimination with complete pivoting may 
also be used to obtain accurate rank- revealing decompositions of graded matrices |17j . 

7.4. Modified one-sided Jacobi algorithm . The heart of the algorithm in ^7, Algorithm 3.1] 
is the modified one-sided Jacobi algorithm, which accurately computes the SVD of matrices of the 
form DB and BD, where D is diagonal and typically highly graded, and B is well-conditioned (see 
|20| . |33) . |24i I25) V Although we focus on the one-sided Jacobi algorithm as applied to G = BD, 
analogous considerations apply to G = DB by replacing G by G*. The one-sided Jacobi algorithm 
works by applying a sequence of Jacobi matrices Ji, . . . , Jm to G from the right (i.e., the same side 
as the scaling, which ensures that components of the right singular vectors are computed with high 
relative accuracy). Each Jacobi matrix J is chosen to orthogonalize two selected columns, and one 
sweep consists of orthogonalizing columns in the order (1, 1), (1, 2), . . . , (1, n), followed by columns 
(2, 3), (2, 4), ... , (2, n), and so on. Sweeps are repeated until all the columns are orthogonal to each 
other to within the bound 

G(Ji---Jm) = VF, ,1/2 ^ne, ii t ^ j. 

\w*Wi\ ' \w*Wi\ ' 

This stopping criterion is used to ensure that even the smallest singular values are computed with 
high relative accuracy. The SVD of G = UTiV* immediately follows by taking Tju = (:,«), 
V = VF/S, and U = (Ji J2 • • • Jm)* ■ 

It will be crucial for the error bounds developed in this paper that the components of the left 
singular vectors of DB (or the right singular vectors of BD) scale in a way similar to D, and are 
computed accurately relative to this scaling. At each step m of the Jacobi algorithm, we write 
( Jq • • • Jyji) G — B„iDjyi, where the columns of have unit P-norm and the matrix is diagonal. 
Defining 

(7.11) vo = max K2 (Bm) , 

l<m<M 

we then have the following result from [33] and |20| . 

Theorem 18. Let G = DB he anxn full-rank, complex-valued matrix, where the diagonal matrix D 
is chosen so that the P-norm of each column of B is unity. Suppose that one-sided Jacobi algorithm 
is used to compute an approximation Ui to the ith left singular vector Ui of G, corresponding to 
singular value 'Su, and the iteration converges after M sweeps. Then the following error hound 
holds on the computed components of Ui : 

Djj ( p{M,n)vl , 2^ 



(7.12) |u. {j) - {j)\ < min ^ / ' % + O {e 

L V ^ii ^jj J \ Tcigap. 

where 

1 ki-o-jl 
relgap, = , _ , 

•tional to M ■ r?!"-^ , a 

value Ttii satisfies 



uj 

p{M,n) is proportional to M ■ r?!"-^ , and vq in defined in il.ll]) . Moreover, the computed singular 



< y^O (e) 
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For convenience, we define 
(7.13) v = p{M,n)vl. 

The following result (see Theorem 3.6 of [34]) will also be needed. 

Theorem 19. Let G = DB, where B is non-singular and D is diagonal. Then the ith left singular 
vector ofG, corresponding to the simple singular value "En, satisfies 

max{||i?-i||,||i?||} . / D,, 



relgapi [ VS~ Djj 

Moreover, there is a left singular vector Ui of the perturbed matrix G+5G = D {B + SB) that satisfies 

max \ \\B^^\f ,k{B)\ r n . . ^/y-) 

\u^U)-u,{j)\ < ^ Ln,m\^,y-^\0{\\6B\\). 

relgap^ Djj J 



8. Appendix: Proof of Proposition [TU] 

8.1. Overview of the proof. Proposition [TO] concerns how graded perturbations of the form 
D {B + 5B) D perturb the singular vectors of DBD, where the diagonal of D > is decreasing 
and B is complex symmetric and non-singular. As in |17) . we analyze the SVD of DBD through the 
LU factorization B — LbL^. In particular, we show that perturbations D {SB) D of DBD result in 
graded perturbations SU of the singular vectors: 

max { jD (SU) E-i/^jj ^ ||^-i (^gjj^ j.i/2|| | < {Lb^Ub) O {\\SB\\) . 

Here csu (Lb) depends on the condition number of Lb- 

We now discuss the basic scaling considerations behind the proof of Proposition [101 First, if 
k{Lb) is not too large, then DBD has the rank-revealing decomposition DBD = XD^Y, where 
X = [dLbD-^) and Y = D^^L^D. Therefore, letting X = QR denote a QR factorization of X 
and setting A = {D-'^RD'^) Y* , it follows that DBD = QD'^A, where 

k{X)<k{Lb), n{Y)<K{LB), ti{A)<K'{LB). 

From theory developed in [34] and [20] (see Lemma [22]) . the matrix of left singular vectors Ui of 
D^A — UiYjU* is graded in the sense that 



(8.1) 



max ■ 



d-^uy}/"^ 



I < max 1 1 



the elements of 



It also can be shown (Lemma [20]) that in the QR factorization QR = {DLbD ^) 
Q are also graded in the sense that 

(8.2) Tiv&y.{\\D-^QD\\ , HDQU^^ || } < h{Lb) ■ 

Finally, from (18. ip and (18. 2p . it is not difficult to show that the matrix U ^ QUi oi left singular 
vectors of DBD is also graded, 

max|||i5'i {QUi) E^'^\ , (g(7/) ^-^/^ || } < /ii {Lb) , 

where the factor /ii (Lb) is defined by 

Ml (Lb) = K(LB)max|||Lg^||^,||L_B||^| . 



CON-EIGENVALUE ALGORITHM FOR OPTIMAL RATIONAL APPROXIMATIONS 



32 



8.2. Preliminary Lemmas. Before proeeeding to the proof of Proposition [TUl we first need several 
lemmas on graded matrices. 

Lemma 20. Let L denote a nonsingular lower triangular matrix, and let D denote a diagonal matrix 
with positive, decreasing diagonal elements. Then if A = DLD~^ has the QR factorization A — QR, 
the elements of Q are bounded by 

m£ix{\\D-'^QD\\ , \\DQD-'^\\} < k{L). 
Proof. Rearranging Qi? = DLD^^, 

D-^QD = L (D-^R-^D) . 

Since R~^ is upper triangular and the elements of D > are decreasing, Lemma [11] implies that 

ll^'^^Q^'ll < \\L\\ < f^iL). 

In the last inequality, we used that ||^^^|| = ||-^-^-^^^|| — ll^ll {s^Ss^in, by Lemma [TT|) . 
To bound ||_D(5_D~^||, we rearrange QR = DLD^^ as 

D-^Q-^D = {D-^RD) L-\ 

Also, since Q^^ — Q*, we calculate that 

\\D-^Q-^D\\ = \\D-^Q*D\\ = [DQD-^] 

= \\{DQD-^)\\. 

Therefore, 

\\{DQD-^)\\ = \\{D-^RD)L-^\\ < \\R\\ < n{L), 

where we again used Lemma [TTl □ 

Lemma 21. Let D and T, denote diagonal matrices, and suppose that the matrices Q and U satisfy 



max{\\D-'^QD\\ ,\\DQD-^\\} < cq, max | 



} < cu. 



Then 



max I L>-i (QC/)Si/2 , D {QU)T.-^/'^ } < cqCj/- 

Proof. We have that 



< 



cqCu- 



The reverse inequality follows from 



< 



CQCu. 



□ 



The following lemma (see Corollary 3.3 and Proposition 3.6 in [34]) shows that the matrix of left 
singular vectors of D^A are graded in a particular way. 

Lemma 22. Suppose that D is a positive- definite diagonal matrix, and A is non-singular. Then if 
D^A has the SVD D^A = UiY^Ur, 



and 



\D-^UiYs\\ < \\A\\ , \\{D^UiJ:-'^)\\ < \\A-^\ 



(C/,),^.|<..(A)min|-^,-^ 
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Lemma 23. Suppose that B is complex symmetric, and has the LU factorization B = LsLg. 
Define 

fii (Lb) = K (LB)max|||ig^||^ , ||Lb||^| . 

Then the right singular vectors of DBD and D {B + SB) D coincide with the right singular vectors 
of D^A and {A + SA), where the matrices A and 5 A satisfy 

\\A\\<\\Lb\\\ \\A-^\<\\L-^'\\\ \\6A\\<^'iLB)0{\\6B\\). 

Proof Note that 

DBD = D{LbLI)D 

= (dlbD-^) d^ (d-^lId) 

= =LiD'^lI. 

Now let Li — QR denote the QR factorization of ii, and define A = (^D^^ RD^) Lj . Then it follows 
that QD^A = LiD^LJ , and the singular values of D^A and DBD coincide. Also, using Lemma [TT] 
we calculate that 

\\A\\<\\Lb\\\ \\A-'\\<\\L-^'\f. 
We now consider D {B + 6B) D. First note that B + 6B has the LU factorization B + SB = 
{Lb + 5Lb){lI+ SUb), where (see [37]) 

(8.3) max{||(5LB||,||,5f/i3||} < \\Lg^\\ k (Lb) \\6B\\ . 
Therefore, 

D{B + dB)D ^ [Li + 5Li) D^ {lJ + dUi) , 

where 

(8.4) ll^Lill = \\DSLbD-^\\ < \\SLb\\ < \\Lg^\\ n {Lb) \\5B\\ . 
In the first inequality above, we used Lemma [TTl Similarly, 

(8.5) ||<5;7i|| < \\Lg^\\n{LB)\\5B\\. 

Now let Li + 5Li have the QR factorization Li + 5Li = Q {R + 6R). Then from [37] . 

\\6R\\ <k(Li)PLi|| < WLg^K^ iLB)\\SB\\. 

Therefore, if we define 

A + SA = {D-^ {R + SR) D'^) {lJ + SUi) 

= A+ (D-^SRD^) lJ + (D-^RD^) SUi, 
we may use (|8.4p and (|8.5p to bound 

< \\Lg'\\K.^LB)\\LB\\\\SB\\ + \\LB\\\\Lg^\\ti{LB)\\5B\\ 
= ^^Lb)0{\\SB\\). 

□ 



Lemma 24. Suppose that B is complex symmetric, and has the LU factorization B 
Define 

fJ-i (Lb) = K (LB)max|||ig^||^ , . 
Then the matrices U and V of left and right singular vectors of DBD satisfy 



^B^B- 



(8.6) 



max I D-if/Ei/^ ^ Df/E-i/z | < (L^) 
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and 

(8.7) 

Also, 



max 



1/2 



DVT. 



-1/2 



} < in [Li 



n . . ' n. 



Proof. Since B is complex symmetric, B has the SVD B = UTV* , V — U. Thus, it suffices to 
consider the matrix U of left singular vectors. 

As in the proof of Lemma we may factor DBD = QD^A, where Q is orthogonal and 

\\A\\<\\Lb\\\ \\A-'\\<\\L-^'\\\ 

Now suppose that D'^A has the SVD D'^A — UiTV* . From Lemma the matrix Ui of left singular 
vectors satisfies 

||i:>"2[/;S|| < ||A|| , \{D^Ui2r^)\ < \\A-^\\ . 



In particular, 



and so 



Z)2 



< \\A-'\ 



(C/0,,| < max{||A-i||,P||}min|g,^| 
< max|||L^^|| , llisll^l min "I 



We can write the inequality above as 



max I D-^UiY}/'^ 



} <max{||is'ir,l|iB||'} 



Also, from Lemma the unitary matrix Q from the QR factorization Li = QR satisfies 
(8.8) max{||L>"igi:>|| , ||L>gL>"i||} < k{Lb) 

Since DBD = QD^A, Lemma [H] shows that the left singular vectors U ~ QUi of DBD satisfy 



(8.9) 



DUT-^I"^ 



} <Mi {Lb): 



where /ii {Lb) is defined in the statement of the lemma. 



□ 



Lemma 25. Let B he complex symmetric, with an LU factorization B — LbL]^. Let D > be a 
positive-definite diagonal matrix, and define 

DiBiDi = 

Then the ith eigenvector of DiBiDi, corresponding to eigenvalue Tu, may be chosen so that the 
following component-wise bounds hold: 


















I) 


V B* 








I) 



(8.10) \x,{j)\<V2fii{LB)mh 
The following norm-wise bound also holds: 



D 



n 



DiXiK, 



-1/2 



> IISII 



-1/2 
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Proof. Since B is complex symmetric, B has the SVD B ~ U^V* , V — U. Therefore, DiBiDi has 
the eigenvalue decomposition 



DiBiDi 



J_fU U 
V2\U -U 



S 
-E 



1 



U U 



V2\U -U 



From (|8.6p in Lemma [Ml 

(8.11) max|||L)-iC/Si/2 , DL/S^^/^ | < k (is) max | f , ||Lsf } 

Since 



lEl^/^ 



,1/2 



1 / DC/ISI^/^ L>C/|S|^/^ 



\/2 V Dt/lS 



1/2 



1/2 



we have from (|8.1ip that 

\d^^Ui |Si|'/'|| < ^/2K(iB)max{||L5lf , • 
An analogous calculation shows that 

|l?i[/i|Eir'/'|| < V2«;(Ls)max{||i^if 

The previous two inequalities show that the ith eigenvector Xi of DiBiDi (e.g., the ith column of 
Ui) satisfies the component- wise bound 



D 



33 



We now consider IjDiSiH /y/Y^u. Since a;,; = (l/\/2) [uiUi]^ , we have from the norm-wise bounds 
in Proposition [TOl that 



-1/2 



1 



n — -^-1/2 



1/2 



DiUiK 



-1/2 



> IISI 



-1/2 



□ 



8.3. Proof of the main result. 



Proof. From Lemma [221 the right singular vectors of DBD and D {B + SB) D agree with the right 
singular vectors of D^A and {A + SA), where the matrices A and E satisfy 

(8.12) \\A\\<\\Lb\\\ < , \\6A\\<>^^iLB)0{\\SB\\). 

Also, Lemmas [52] and (Ml ensure that the ith left and right singular vectors Ui and Vi of D^A are 
bounded by 



D 



33 



• 13) max{|ui (j)| , \vi < fii {Lb) min ■ 

L V ^33 

To prove the component-wise bounds, we now proceed as in Theorem 2.21 of [20] (except now we 
may use the above component-wise bounds for both Vi (j) and Ui (j)). Namely, let SA — r/E, where 
rj = \\SA\\, and define 



Z?2 
/ 



El 



E* 
E 



Ai 



A* 
A 



1 



Vi 



V2 V 
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Let xf (77) denote the eigenvectors of G {rj), 
G{v) = 

Then by standard perturbation theory, 



{A + riE)*D'^ 

D^{A + riE) 



Di {Ai+7]Ei)Di. 



and so 



Now, 



Xa ) D-\E-\ D^Xi. I , o\ 



±k^±i 



±k^±i 



--yD^E^D^xi\ < ^\{v,,D^E*Uk)' ^ 



,ED\k)\ 



= ^\{D\„E*Uk)\ + l\{u„ED\k)\ 



< 



D^vk 



Now, we have that D'^Avk — UkUk, and so 

<Jk - \\D^Avk\\ > \\A-^f^ \\D^Vk\ 
Similarly, ||Z?^?;J| < (Ti 11 A^-'- 11 . Therefore, 



\D\k\\ < (Jk \\A-^ 



(8.14) 



[xfY'DiEiDix 



Finally, from ((51^ and ((gT^ . 



< 



(xf)* DiEiD^xf 



< 0(77)||A-i||Aii(£B) E 



< 0(?7)||L^M| ^1 (Lb) (max 



2 ,^yf (^i+<7k 



D 



n 



/Y^i Djj 



\ k^i |(Tj — fjfcl 



33 



< Oi\\SB\\) \\LgH\\i{LB)K^LB) (max 



D 



33 



where we used (|8.14p in the second to last inequality and (I8.12p in the last inequality. 

We now prove the norm-wise bound on Dui/Yin. To do so, note that, since B is complex sym- 
metric, the SVD of B may be written as _B = UY^U"^ (e.g., B has a Takagi factorization). Now 
suppose that DBDul = UiUi. Then the bound follows from 

(J, = \{DBDul,Ui)\ = \{BDlM,Dur)\ < \\B\\ \\Du.,f . 

Finally, the bounds on the ratio Ya/Yia of the singular values of DBD and D {B + 6B) D follow 
from Theorem 4.1 of [T7]. □ 
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